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Preface 


The  purpose  of  the  study  was  to  investigate  the  ability  of  a  Kalman  filter 
to  determine  the  relative  position  of  satellites  operating  in  a  cluster  environment. 
The  investigation  is  a  follow-on  to  a  study  presented  by  Captain  Michael  Ward  in 
December  1988  and  includes  the  addition  of  other  forms  of  measurement  data  to  the 
filter  and  a  study  of  the  observability  of  the  state  components. 
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Abstract 

The  relative  position  determination  of  a  cluster  of  satellites  in  a  near  circular 
orbit  was  investigated  in  a  thesis  by  Captain  Michael  L.  P.  Ward  in  December  1988. 
His  investigation  involved  the  use  of  dynamics  based  on  the  Clohessy-Wiltshire  equa¬ 
tions  and  an  on-board  estimator  based  on  the  U-D  covariance  factorization  version 
of  the  Kalman  filter.  The  initial  performance  results  proved  favorable,  with  the  filter 
meeting  the  required  25  meter  accuracy  in  all  test  cases.  The  purpose  of  this  thesis  is 
to  validate  the  test  results  achieved  by  Captain  Ward  and  to  investigate  the  ability  of 
the  filter  to  determine  the  relative  position  of  the  satellites  to  a  higher  degree  when 
the  filter  in  use  is  not  resident  on  that  satellite.  This  investigation  included  the  use  of 
additional  measurement  data  from  other  satellites  in  the  cluster  to  the  Kalman  filter 
for  processing  in  the  update  cycle.  When  it  was  determined  that  the  measurements 
were  not  allowing  satellite  1  to  update,  the  observability  of  the  components  of  the 
state  vector  were  computed  and  the  results  discussed. 
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INVESTIGATION  OF  THE  OBSERVABILITY 
OF  A  SATELLITE  CLUSTER 
IN  A  NEAR  CIRCULAR  ORBIT 


I.  Introduction 

The  construction  of  a  group  of  satellites  into  clusters  has  been  under  investi¬ 
gation  in  the  area  of  satellite  communications  for  severed  years  (5:831).  A  cluster 
is  defined  as  an  array  of  two  or  more  satellites  located  in  a  specified  orbital  area, 
that  appear  and  operate  as  one  to  an  Earth-based  user.  The  amount  of  separation 
allowed  between  satellites  is  determined  by  its  application  and  the  accuracy  of  the 
orbital  control  method.  The  benefits  to  employing  this  strategy  are  numerous.  By 
dividing  the  work  load  between  satellites,  the  entire  system  can  provide  a  service 
greater  than  each  individual  satellite.  The  entire  system  is  more  reliable  due  to  the 
availability  of  in-orbit  spares  and  survivability  is  enhanced  because  destroying  one 
satellite  does  not  destroy  the  entire  system. 

There  are  several  problems  that  need  to  be  addressed  in  order  for  this  cluster 
strategy  to  work.  These  areas  include  cooperative  orbit  control  strategies,  cluster 
geometry  and  orbit  determination.  This  thesis  is  concerned  only  with  the  problem  of 
accurately  determining  the  relative  position  of  each  satellite.  The  stochastic  nature 
of  this  position  determination  problem  has  previously  been  address  by  author  John 
Murdoch  in  several  articles  (7:1001-1018).  His  method  of  analysis  used  a  least  squares 
algorithm,  which  is  less  appropriate  than  a  recursive  filter  for  on-board  applications, 
and  his  cluster  size  was  limited  to  two  satellites.  This  thesis  determines  the  position 
of  up  to  ten  satellites  in  a  cluster  and  uses  a  recursive  Kalman  filter  to  estimate  the 
position  of  the  satellites  (5:1006).  The  approach  presented  in  this  study  takes  the 
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pulses  used  by  each  satellite  to  synchronize  their  on-board  clocks  to  provide  input  to 
a  Kalman  filter  for  satellite  relative  position  determination.  This  filter  was  designed 
for  satellite  1  but  can  easily  be  adjusted  for  use  with  any  of  the  members  of  the 
cluster  by  having  each  satellite  think  it  is  number  1.  The  assignment  of  numbers 
to  the  satellites  is  an  arbitrary  designation  since  the  cluster  geometry  is  randomly 
determined  by  Equation  (3). 

The  proposed  use  of  this  Kalman  filter  is  in  a  cluster  designed  as  a  space  based 
radar.  In  order  for  the  radar  to  form  a  cohesive  image,  the  relative  position  of  each 
satellite  needs  to  be  known  to  at  least  one  quarter  of  a  wavelength.  At  the  one  quarter 
limit  however,  the  resultant  fogginess  and  loss  of  contrast  are  appreciable  and  can 
affect  the  visibility  of  very  fine  detail  (2:443-444).  Since  the  military  requirements 
of  a  phased  array  radar  may  include  the  need  for  fine  detail,  the  tolerance  will  be 
decreased  to  one  tenth  of  a  wavelength  or  21.6  meters.  The  method  used  to  obtain 
this  value  is  presented  in  Appendix  A. 

The  study  presented  by  Captain  Michael  Ward  proved  that  the  concept  of 
using  a  recursive  filter  to  determine  relative  position  is,  under  initial  testing,  viable 
(8.T-1-4-1).  The  purpose  of  this  thesis  is  to  continue  the  testing  of  the  U-D  covariance 
factorization  Kalman  filter  (3:392-400),  to  include  determining  filter  error  statistics 
and  to  perform  an  investigation  of  the  discrepancy  between  the  results  achieved 
for  satellite  1  and  those  achieved  for  other  members  of  the  cluster.  The  measure 
of  performance  will  be  the  comparison  of  the  average  of  the  standard  deviations 
extracted  from  the  eigenvalues  of  the  position  covariance  matrix  to  the  statistics 
of  the  true  error.  The  true  error  is  defined  as  the  root-mean-square  of  the  relative 
position  components  computed  from  the  estimated  state  provided  by  the  filter  and 
the  actual  state  provided  by  the  truth  model.  For  this  performance  investigation  the 
cluster  was  placed  in  near-circular,  low-earth  orbit.  (8:1-2). 
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II.  Background 


Since  this  thesis  is  a  continuation  of  Captain  Michael  Ward’s  work,  the  de¬ 
velopment  of  the  initial  Kalman  filter  and  truth  model  will  be  represented  in  less 
detail.  For  a  more  detailed  analysis  of  the  filter  development,  the  reader  should 
consult  the  Ward  thesis  (8).  The  truth  model  was  based  on  an  array  of  from  one  to 
ten  satellites  orbiting  the  earth  with  the  same  orbital  period.  The  truth  model  was 
limited  to  the  two  body  equations  of  motion  and  had  as  output  the  true  state  xf(t), 
whose  components  are  the  relative  position  and  velocity  of  each  satellite  with  respect 
to  a  rotating  reference  point,  and  the  satellites  relative  position  measurement  zt(t). 
This  relative  measurement  data  was  then  used  by  the  Kiilman  filter  to  determine  an 
estimated  state  x(f).  The  position  elements  of  the  true  state  were  then  compared  to 
those  of  the  estimated  state  for  each  satellite  to  provide  the  true  error  of  the  system 


Figure  1.  Kalman  Filter 

The  true  error  statistics  were  then  compared  to  the  average  value  of  the  stan¬ 
dard  deviations  defined  as  the  (y/ eigenvalues)  of  each  satellite’s  position  covariance 


3 


matrix  to  determine  the  filter’s  performance.  The  covariance  of  the  filter  tells  the 
user  how  well  the  filter  thinks  it  is  performing. 


2.1  Truth  Model 


The  truth  model  provides  the  filter  with  only  relative  information  concerning 
the  position  of  each  of  the  satellites  in  the  cluster.  Since  the  measurement  data 
is  relative,  absolute  position  can  not  be  estimated.  Instead  the  relative  position 
was  determined  with  respect  to  a  rotating  reference  point.  The  earth-centered  frame 
(i,j,  k)  is  considered  inertial  and  the  reference  frame  (x,y,  z )  is  rotating  with  respect 
to  this  inertial  frame  at  a  constant  angular  velocity  u  k.  The  reference  point  is  in 
a  circular  orbit  with  a  radius  of  R  and  has  a  circular  orbital  velocity  of  This 
arrangement  is  shown  in  Figure  2. 

The  position  and  velocity  of  the  reference  point  in  the  inertial  frame  at  the 
initial  time  (t=0)  are: 


rre/ 


-R- 

0 

0 


(1) 


0  1 


Vref 


Vcs 

0 


(2) 


The  positions  of  the  satellites  were  determined  at  random  about  the  reference 
point  but  within  a  specified  radius  (rad)  and  the  velocity  vectors  were  determined 
from  the  Clohessy-Wiltshire  equations  (10:80).  These  equations  have  the  following 
form: 
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Figure  2.  Reference  System 


R ‘ 

n  =  0  +  rii  *  (rad) 

.0. 

(3) 

Vii  =  T](ni-  rPe/  i) 

(4) 

Ui  J  =  (2 /x(i  -  -  v? t  -  uf  fc)’ 

(5) 
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A  A 

Vi  k  =  rjrik 


(6) 


where  i„  the  output  of  a  uniformly  distributed  random  number  generator  with 
outputs  between  —0.5  and  0.5  and  ij  is  defined  as  the  mean  motion. 

The  input  for  the  Kalman  filter  can  be  formed  using  this  position  and  velocity 
information.  The  true  initial  state  xt(t  =  0)  is  defined  as  the  position  and  velocity  of 
each  satellite  in  the  rotating  frame.  Since  the  coordinate  frames  are  initially  aligned, 
the  position  of  each  satellite  in  the  rotating  frame  is: 

ri[rot]  =  (ri  —  rref)jix  (7) 

The  velocity  expressed  in  term  of  relative  position  is: 

^[rot]  =  (vi  ~  Vref)yix]  —  UJ  X  7\[rot]  (8) 

Ward  discovered  in  his  thesis  that  the  absolute  downrange  position  y  was 
unobservable.  The  downrange  component  was  removed  from  satellite  1  and  the 
state  of  satellites  2  through  s  were  modified  to  include  the  relative  (Ay,)  rather  than 
the  absolute  downrange  component,  where  Ay*  is  defined  as: 

A  yi  =  yi  -  yi  (9) 

with  i  =  2,3,  ...,s.  Figures  3  and  4  show  the  x  and  y  components  of  the  state  with 
and  without  this  modification.  The  z  component  is  normal  to  the  plane  of  the  page 
and  was  not  affected  by  the  modification. 

With  this  modification  the  state  vectors  in  the  rotating  reference  frame  are  defined 
as: 
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Figure  3.  The  Original  Component  Representation  of  the  Cluster 


Figure  4.  The  Modified  Component  Representation  of  the  Cluster 
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The  relative  measurement  vector  has  the  following  form: 


'  I  n  -  r3  | ' 

Zt  =  •  +  Ut  (11) 

-lri  |. 

where  ut  is  a  zero-mean,  white  Gaussian  noise  diagonal  matrix  with  a  covariance  of 
Rt  (3:330).  The  values  for  the  covariance  Rt  were  determined  by  the  Ward  thesis 
to  be  .01  meters.  This  value  was  probably  obtained  from  the  satellite  with  largest 
true  error  calculated  once  the  filter  was  tuned.  Statistically,  this  allows  the  adver¬ 
tised  error  in  the  measurements  to  match  the  actuals  error  in  the  position  of  the 
satellites.  The  noise  ut,  represents  the  errors  in  computing  the  range  measurements 
from  one  satellite  to  the  others.  The  major  errors  result  from  clock  synchronization 
differences  between  satellites  as  well  as  resolution  limitations  of  the  on-board  clock 
(7:1016).  These  errors  were  also  considered  to  be  independent  from  measurement 
to  measurement  and  were  therefore  grouped  together  as  one  error  term  along  the 
diagonals  of  the  measurement  covariance  matrix. 

The  solution  to  the  Kepler  problem  is  used  to  determine  the  future  position  and 
velocity  of  each  satellite.  The  Kepler  problem  is  solved  using  the  /  and  g  equations 
in  terms  of  the  eccentric  anomaly  E  (1:219).  These  equations  have  the  following 
form: 
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(12) 


/  =  1 - (1  —  cos  A E) 

ro 

=  t  ~  ]fj(AE-*inAE) 

y/jlaamAE 


f  = 


rr0 


3  =  1 - (1  —  cos  A E) 


where  A E  and  r  are  defined  as: 


(13) 

(14) 

(15) 


A  E  =  Ef  -  E0 

(16) 

=  a(l  —  ecos  Ef) 

(17) 

The  value  for  Ef  was  determined  through  the  use  of  a  Newton  iteration  method 
since  the  eccentric  anomaly  does  not  have  a  closed  form  solution  (9:61).  A  more 
detailed  description  of  this  process  is  presented  in  the  Ward  thesis  (8:2-6-2-7).  Once 
the  f,g,f  and  g  variables  have  been  computed  the  new  position  and  velocity  can  be 
determined  using  the  following  equations: 


r(t)  =  f  ro  +  g  vo 

(18) 

v(t)  =  fr0  +  gv  o 

(19) 

Once  the  position  and  velocity  of  each  satellite  in  the  inertial  frame  have  been 
determined,  the  relative  position  to  the  reference  point  can  be  found.  The  first  step 
is  to  determine  the  inertial  position  and  velocity  of  the  reference  point.  Defining  9 
as  u>t,  the  position  and  velocity  vectors  are: 
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’  R  cos  0 " 

rT  e/[/i*](0  =  -ft  sin  0  (20) 

0 

— vet  sin0' 

VTef[fix](t)  =  uCJcos0  (21) 

0 

The  next  step  is  to  determine  the  relative  distance  between  the  reference  point 
and  each  satellite  in  the  inertial  frame.  Using  the  position  and  velocity  of  the  satel¬ 
lites,  computed  from  the  eccentric  anomaly ,  and  the  position  and  velocity  of  the 
reference  point  listed  above,  the  relative  distance  becomes 

rrel[fix)  =  K*)  -  rref(t)]  =  Tj  i  +  r2  j  +  r3  k  (22) 


Vrellfiz]  =  [v(0  -  Vref(t)]  =  *  +  V2  j  +  V3  k  (23) 

Since  the  inertial  and  rotating  reference  frames  are  no  longer  aligned,  a  co¬ 
ordinates  transformation  must  be  performed  to  obtain  the  relative  distance  in  the 
noninertial  frame.  Performing  this  transformation  yields: 

ri  cos  9  -f  r2  sin  6  1  x 

—r\  sin  9  +  r2  cos  9  y 


rrei[rot] 

vr*J[rot] 


cos  9  +  tij  sin  0  +  u>  rj 
— vj  sin  0  +  v 2  cos  0  —  u>  rj 
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2.2  Kalman  Filter 


The  truth  model  presented  in  the  last  section  provides  the  actual  state.  The 
purpose  of  the  Kalman  filter  is  to  estimate  this  state  to  a  predetermined  degree  of 
accuracy.  This  filter  performs  two  functions:  the  first  is  to  propagate  the  estimated 
state  and  its  covariance  forward  in  time,  and  the  second  is  to  update  the  state  and 
covariance  using  the  measurement  data  provided  by  the  truth  model.  This  cluster 
system  is  represented  in  the  form  of  a  linear  stochastic  differential  equation  which 
describes  the  state  propagation,  with  discrete- time,  noise  corrupted,  nonlinear  inputs 
as  available  data  (8:2-10). 

2.2.1  State  Propagation  and  Update  . 

Since  the  dynamics  used  to  propagate  the  state  are  discrete-time  and  linear, 
the  underlying  dynamics  model  has  the  following  form: 

x(f.)  =  i)  +  (25) 

where  wj  is  a  discrete-time,  zero-mean,  white  Gaussian  noise  sequence  with  covari¬ 
ance  equal  to  the  dynamics  driving  noise  Qj  at  time  equals  U,  and  Gj  is  assumed  to 
be  the  identity  matrix  for  an  equivalent  discrete-time  model  (3:220).  The  $  matrix 
represents  the  state  transition  matrix  of  the  system  and  will  be  defined  in  more  detail 
later  in  this  section.  This  equation  leads  to  the  standard  filter  propagation  equation 
for  the  Kalman  filter  defined  as  (3:220): 

x(ti  )  =  $(t»,  *«-i)*(f£.i)  (26) 

where  ( t~ )  and  (t+)  represent  the  state  before  and  after  the  update.  For  convenience, 
they  will  be  represented  as  (  — )  and  (-f)  for  the  remainder  of  this  thesis.  Each 
satellites  $  matrix  is  determined  from  the  Clohessy- Wiltshire  equations  of  motion, 
with  7)  representing  the  mean  motion  (10:80).  These  equations  axe: 


11 


x  —  2rjy  —  3rj2  x  =  0 


(27) 


y  +  2rj  x  =  0  (28) 

z  +  rj7  z  =  0  (29) 

The  above  equations  are  valid  for  small  displacements  in  the  radial  and  out-of-plane 
directions  but  remain  correct  for  large  changes  in  the  in-track  direction  (10:80). 
These  equations  of  motion  are  integrated  about  the  initial  conditions  xo,x0,yo,y0, 
etc.  at  t  =  0  to  obtain  the  following  position  and  velocity  solutions  (8:2-12): 

(2  \  xo  2 

-  j/o  +  3  x0  cos  rjt-\-  —  sin  rjt  +  4  x0  4-  -y0  (30) 

V  J  V  V 


y(0  =  yo  -  (3  2/0  +  6tj  x0)t  4-  (  ^2.  +  6x0  I  sin  rjt  +  ^^cos  r/t  -  (31) 

\  V  J  V  V 

z(i)  —  zq  cos  rjt  +  —  sinr;<  (32) 

V 

x(t)  =  (2  yo  +  3  rj  x0)  sin  rjt  +  io  cos  rjt  (33) 

y(t)  =  -3  yo  —  6  rj  xo  +  (6  rj  xo  +  4  y0)  cos  rjt  -  2  x0  sin  rjt  (34) 

z(t)  =  — zqtj  sin  rjt  +  Zq  cos  rjt  (35) 
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Using  these  equations,  the  following  $  matrix  can  be  formed  for  a  satellite 
with  six  states  represented  by  x,y,z,x, y,  and  z : 


4  —  3  cos  ip 
6(sinV'  —  ip) 


*i  = 


0 

3  77  sin  ip 

6  Tj  (cost p  —  1) 
0 


0 

1 

0 

0 

0 

0 


0 

0 

cos  ip 
0 
0 

—77  sin^ 


«in  t/i 

1 

5(1  -  cos  VO 

0 

ij  (cos  ip  1) 

4  sin 

0 

0 

0 

tin 

*J 

cos  ip 

2  sin  ip 

0 

—2  sin  ip 

—3  +  4  cosip 

0 

0 

0 

cos  ip 

(36) 


where  ip  =  rjt.  Since  the  state  vectors  were  modified  to  include  Ay;  instead  of  y;, 
Equation  (31)  was  substituted  into  Equation  (9).  For  t  =  2,3,  ..,s  this  modification 
yields: 


Ay.(+) 


(6  (sin  ip  -ip))  (x  1-  Xi)  +  Ay;  + 

(^sin^- y)  fa-Vi) 


-  (cosip  —  1) 
V 


(l!  -  Xi)  + 


(37) 


where  the  components  on  the  right  hand  side  of  the  equation  are  evaluated  at  the 
previous  sample  time.  Each  satellite’s  $  matrix  can  now  be  formed.  $1  is  reduced 
from  a  6-by6  to  a  5-by-5  by  matrix  with  the  row  and  column  corresponding  to  the  y 
component  removed  (8:2-19): 
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$1  = 


4  —  3  cos  V* 

0 

3  rj  sin  ip 
677  (cos  ip  —  1) 


cos  ip 


—77  sinip 


**  ;(!-«.«  o  - 

0  0  ^ 

*1 

cos  ip  2  sin  ip  0 

—2  sin  ip  — 3  +  4cosV>  0 

0  0  cos  *0 . 


The  $  matrices  for  the  remaining  satellites  we: 


4  —  3  cos  ip  0  0 

—6  (sin  ip  —  ip)  1  0 

0  0  cos  ip 

3  rj  sin  ip  0  0 

677  (cos  ip  —  1)0  0 

0  0  —77  sin^» 


|(1 -co.*)  0 

(cosip  -  1)  — -sin^»  +  ^  0 

ti  '  '  '  n  1  fi 


cos  ip 
-2  sin  ip 


2  sin  ip 
—3  +  4  cos  ip 


cos  ip  _ 


where  the  sign  of  components  in  the  second  row  have  been  changed  except  for  the  y 
component  due  to  the  way  Ay*  was  defined. 

The  overall  $  matrix  for  this  system  has  non- diagonal  terms  due  to  the  inter¬ 
action  between  the  y  coordinates  of  each  satellite.  The  off-diagonal  elements  are: 

0  0  0  0  O' 

6 (sin ip  — ip)  0  jj  (cos ip  —  l)  ;j-sin i>  —  ^  0 

0  0  0  0  0 

in  =  (40) 


The  overall  system  has  the  following  form: 
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'$1  0  ...  0  " 

’®l(+) 

$ai  $a  •  •  •  0 

•  •  •  • 

*2(  +  ) 

•  •  •  « 

•  «  . 

0  0 

**(+) 


(41) 


The  state  is  propagated  from  t  =  0  to  the  desired  St  before  any  updates  are 
made,  therefore  the  initial  value  for  the  state  must  be  determined.  The  initial  value 
for  the  estimated  state  x0(f  =  0)  was  set  equal  to  the  initial  value  from  the  truth 
model  xt(t  =  0)  with  the  position  measured  in  kilometers  and  the  velocity  measured 
in  kilometers/second. 


Once  the  state  has  been  propagated  through  the  first  time  interval,  it  must  be 
updated.  Since  the  measurements  are  discrete-time  and  nonlinear,  the  generic  model 
for  the  measurements  has  the  following  form: 


z(U)  =  A(x(ti),  U]  +  v(U )  (42) 

where  h  is  the  filter’s  nonlinear  estimate  of  the  measurements  and  v  is  a  white 
Gaussian  noise  sequence  of  zero  mean  and  covariance  of  Rt  which  was  described 
previously  (4:40).  This  leads  to  the  extended  form  of  the  Kalman  filter  update 
equation  (4:44): 


x(+)  =  i(— )  +  *■{*  —  >>  [*(—)]}  (43) 

where  K  is  the  gain  which  will  be  determined  once  the  covariance  has  been  updated, 
and  h[x(— )]  is  the  filter’s  nonlinear  estimate  of  the  satellite’s  relative  range  before 
the  state  has  been  updated.  The  estimate  of  relative  range  has  the  following  form: 
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(44) 


'  x2f  +  Ayf  +  (zi  -  z2f ' 

= 

•  hx  ■ 

.  ^/(ii  -  x,)2  +  A y)  +  (zi  -  z,)2 . 

2.2.2  Covariance  Propagation  and  Update  . 

In  order  to  obtain  numerical  stability  with  single  precision,  the  U-D  covariance 
factorization  version  of  the  Kalman  filter  was  used  to  propagate  and  update  the 
covariance  (3:392-399).  This  form  of  the  filter  was  necessary  to  guarantee  the  pos¬ 
itive  semi-definiteness  of  the  computed  covariance  matrix  and  to  achieve  numerical 
stability  (8:2-15). 

Double  precision  accuracy  is  obtained  with  single  precision  by  factoring  the 
covariance  into  two  forms,  the  propagation  of  the  covariance  represented  by  P(— ) 
and  the  updating  of  the  covariance  represented  by  P(+): 

/>(-)  =  £/(-)£(-)  VT(-)  (45) 

P(+)  =  U(+)D(+)U*(+)  (46) 

where  the  U  matrix  is  upper  triangular  and  unitary  and  the  D  matrix  is  diagonal. 

The  first  step  in  the  algorithm  is  to  determine  the  initial  values  for  the  two 
matrices  U  and  D  by  using  the  existing  initial  covariance  Po,  which  is  a  diagonal 
matrix.  This  initial  covariance  contains  the  steady  state  values  determined  from  the 
initial  testing  accomplished  by  Captain  Ward,  which  were  1  x  10-6  km  for  position 
and  1  x  10-ia  ^  for  velocity.  This  initial  factorization  is  obtained  using  the  following 
equations: 

Dnn  =  Pnn  (47) 


16 


[1  i  =  n 

Uin  =  {  „  ,  „  (48) 

1  Pin/ Dnn  t  =  n  -  1,  n  -  2,  . . . ,  1 

where  n  represent  the  total  number  of  states  in  the  system  (6  *  s  —  1).  As  described 
previously,  each  satellite  has  six  states  except  for  satellite  1  which  has  five.  The 
remaining  columns  in  the  matrices  are  computed  as  follows: 


Dii  =  Pa  -  E  0“  V* 

k=j+l 


U*  =  1 


Pij  Z)k=j-i  Dkk  U\k  / D ii  i  —  j  —  1)  j  —  2, . . . ,  1 


Now  that  the  initial  values  have  been  determined,  the  covariance  can  be  prop¬ 
agated  to  the  first  update  time.  Using  the  process  above,  the  value  for  U(+)  =  U0 
and  D(+)  =  D0  have  been  determined  and  are  used  in  forming  two  additional 
matrices.  The  first  matrix  is  n-by-2n  and  is  designated  as  Y(-): 

r(-)  =  [$U(+)  |  Gd]  (51) 

where  $  has  been  previously  defined  and  Gd  again  is  an  identity  matrix  because 
the  model  is  an  equivalent  discrete-time  representation  of  a  continuous-time  system 
(3:337).  The  second  matrix  is  a  2n-by-2n  and  is  designated  D{— ): 


D(-)  = 


D(+)  0 

0  Qd 


The  propagation  begins  by  transposing  the  Y{— )  matrix 
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yT(~)  =  \a2a2  ...  an] 


(53) 


and  then  calculating  the  following  relationships  for  k  =  n ,  n  —  1 , . . .  ,  1 


ck  —  D(  —  )ak  (cjk  —  Ajj(— ) ai*i  j  -  l,2,...,2n) 

Dkk(-)  =  al  Cfc 

=  ck/Dkk(-)  (54) 

£fyk(-)  =  aj  4  j  =  1,2,...  fc  —  1 

a,j  *  Oj  Uj  it(  )  flfc  j  —  1,2, . . . ,  fc  1 


Now  all  the  information  has  been  calculated  to  propagate  the  covariance  to  the 
new  time  using  Equation  (45).  The  next  step  is  to  update  the  covariance.  Several 
terms  need  to  be  defined  for  these  calculations.  Let  the  linearized  measurement 
vector  Hi  be  defined  as: 


H-Sh\ 

H~di 


(55) 


Using  this  definition  the  overall  system  H  matrix  has  the  following  form: 


'  Hx  Ht  0  ...  O' 

H2  0  H2  ...  0 


H.- ,  0  0  ...  H »_i 


(56) 


where  the  individual  matrices  for  i  —  1, 2, . . . ,  s  —  1  are  defined  as  follows: 


=  0  0  0 


(57) 
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Hi  =  HtT11)  ^  0  0  0]  (58) 

The  H  matrix  is  divided  into  individual  row  vectors  and  these  are  designated 
as  Ht.  The  variable  R\  is  defined  as  the  individual  diagonal  elements  in  the  data 
covariance  matrix  Rf,  which  was  assumed  to  be  equal  to  Rt  for  this  thesis.  Now  that 
all  the  terms  have  been  defined  the  update  can  be  computed.  For  z  =  1, 2, . . . ,  s  —  1 

/  = 

vi  ~  Djj(-)  fj  j  —  1,2,...  ,n  (59) 

a0  =  Ri 

Then  for  k  =  1,2, . . .  ,n 

<*k  —  CLk- 1  +  fk  Vfc 

Dkk(+)  =  Dkk(—)ak~i/ak 

bk  *-  vk 

Pk  =  -fk/dk- 1  (60) 

Ujk(+)  =  Ujk{—)  +  bjpk  j  =  1, 2,  ...,(&  —  1) 
bj  «-  bj  +  Ujk(—)vk  j  =  1,2,  1) 


Using  the  values  for  £/(+)  and  Z?(-f)  computed  above  the  covariance  can  be 
updated  using  Equation  (46).  The  gain  of  the  filter  is  defined  as: 

On 

with  b  composed  of  the  most  resent  components  calculated  in  Equation  (60).  The 
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state  update  can  now  be  determined  using  Equation  (43).  This  entire  process  of 
propagating  and  updating  the  state  and  covariance  is  repeated  at  a  time  interval 
specified  by  the  user  until  a  final  time  is  reached. 

In  summary,  this  chapter  developed  the  truth  model  based  on  the  two-body 
equations  of  motion.  This  model  consisted  of  the  relative  position  and  velocity  com¬ 
ponents  of  each  satellite  with  respect  to  a  rotating  reference  point.  A  U-D  covariance 
factorization  version  of  the  Kalman  filter  was  used  to  determine  the  estimated  po¬ 
sition  and  velocity  of  each  satellite  using  the  measurement  data  provided  by  the 
truth  model.  A  modification  to  the  y  component  of  both  the  truth  model  and  filter 
became  necessary  due  to  observability  problems.  The  modification  included  replac¬ 
ing  the  absolute  y ;  component  of  satellite  2  —  3,  with  a  relative  Ay;  component  by 
removing  the  y  component  from  satellite  1.  The  state  for  satellite  1  was  reduced 
to  five  components  with  the  other  members  of  the  cluster  retaining  six  components. 
The  next  step  will  be  to  test  the  filter. 


20 


III.  Initial  Testing 


Now  that  the  true  state  and  estimated  state  have  been  computed,  the  true 
error  et  of  the  position  components  can  be  obtained.  The  errors  in  the  velocity  will 
not  be  computed  since  the  purpose  of  this  thesis  is  to  determine  the  accuracy  of  the 
filter’s  estimate  of  the  relative  position  of  each  satellite  in  the  cluster.  The  error  in 
the  velocity  will  cause  the  filter’s  computed  cluster  to  separate  and  this  separation 
will  appear  as  errors  in  the  position  components.  Since  the  elements  of  the  true 
state  and  estimated  state  are  of  the  same  type,  for  t  =  1, 2, 3, . . . ,  s  the  true  error  is 
defined  as: 


et  =  [( Xi  -  xtif  +  ( yi  -  yti)2  +  (i<  -  *ti)2]  *  (62) 

where  the  (yi  —  yti)2  component  is  missing  for  satellite  1.  Once  the  true  error  has  been 
computed  it  can  be  compared  to  the  estimated  error  which  is  defined  as  the  average 
standard  deviation  of  the  position  components  computed  from  the  covariance,  which 
again  is  a  measure  of  how  well  the  filter  thinks  it  is  performing.  The  estimated  error 
is  computed  using  the  following  equation: 


(Tavg  =  +  O ]  (63) 

where  ul  y<J  are  the  eigenvalues  of  the  position  covariance  matrix.  For  Captain 
Ward’s  thesis  the  true  error  was  compared  to  the  maximum  standard  deviation 
(ffmax)  of  the  covariance,  which  was  defined  as  the  square  root  of  the  largest  eigen¬ 
value  of  the  position  covariance  matrix.  The  comparison  was  changed  initially  in  a 
attempt  to  determine  if  the  curves  for  the  true  error  could  be  made  to  match  the  av¬ 
erage  values  of  the  covariance  position  components  to  a  higher  degree.  This  assumed 
that  the  filter  was  performing  properly,  which  was  later  discovered  to  be  untrue.  Any 
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further  investigations  into  this  problem  should  revert  back  to  a  comparison  of  the 
maximum  standard  deviation. 

S.l  Sample  Runs 

Captain  Ward  tested  the  filter’s  accuracy  by  accomplishing  single  performance 
runs,  which  involve  using  only  one  seed  value  in  the  random  number  generator. 
During  the  runs,  the  number  of  satellites  in  the  cluster  was  varied  and  different 
values  of  dynamics  driving  noise  covariance  Qj  were  used  to  tune  the  filter.  The 
filter’s  relative  measurement  covariance  Rf  was  equal  to  Rt  for  all  time  and  was 
based  on  an  assumed  range  error  of  .01  meters.  The  cluster  radius  was  500  meters 
at  a  reference  altitude  of  1000  kilometers  above  the  earth  with  an  orbited  period  of 
approximately  6300  seconds  (8:3-2-3-3). 

Initially  the  filter’s  stability  was  verified  by  Captain  Ward,  and  runs  were 
conducted  to  determine  if  the  true  error  was  within  the  required  accuracy.  The  filter 
performed  within  the  desired  accuracy  in  all  cases,  including  the  addition  of  noisy 
data  to  the  estimator.  The  maximum  transient  response  of  the  estimated  error  was 
periodic  with  a  peak  at  every  half  of  an  orbit  of  approximately  2  meters  and  the 
maximum  value  of  the  true  error,  which  was  also  periodic,  was  approximately  .01 
meters  for  satellite  1  and  close  to  zero  for  the  other  satellites.  When  noise  was  added 
to  the  data,  the  true  error  as  expected,  increased  to  around  .5  meters  and  in  the 
case  of  satellite  1  this  error  grew  to  exceed  the  filter’s  computed  error  statistics  (tr). 
Through  the  addition  of  more  dynamics  noise  in  the  tuning  of  the  filter,  the  true  error 
decreased  to  a  value  below  the  filter’s  estimated  error.  The  first  task  of  this  thesis  was 
to  verify  the  results  obtained  by  Captain  Ward.  The  code  was  obtained  and  it  was 
determined  that  it  had  been  developed  and  run  on  a  Digital  Equipment  Corporation 
VAX-11-785  with  a  VMS  operating  system.  The  code  ran  without  incident,  but 
the  results  obtained  were  not  the  same  as  the  Ward  thesis.  The  results  of  satellite 
1  showed  the  largest  deviation  with  large  peaks  in  the  true  error  occurring  in  the 
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initial  test  run  which  contained  only  two  satellites  and  no  driving  noise  covariance 
Qd ■  Since  the  Ward  thesis  did  not  specify  an  update  time,  the  code  was  rerun  with 
updates  every  1  minute.  The  results  still  did  not  change.  At  this  point  it  was 
assumed  this  was  not  the  latest  version  of  the  code.  Since  the  state  of  satellite  1  was 
decreased  to  five  states  in  the  last  version  of  the  filter,  this  investigation  searched 
for  errors  in  the  computation  of  the  estimated  state  for  satellite  1.  This  proved  to 
be  the  problem.  The  computation  of  the  true  error  for  satellite  1  was  accounting 
for  three  position  components,  not  two.  Several  small  errors  were  also  found  in  the 
computation  of  the  initial  values  for  the  Uq  and  D0  matrices  and  several  matrices 
whose  values  were  computed  in  the  filter  were  initialized  to  zero  to  insure  their 
values.  Later  in  the  testing  it  was  discovered  that  during  the  measurement  update, 
the  components  of  the  {/(— )  matrix  were  being  updated  with  each  measurement, 
but  components  of  the  U(—)T  matrix  were  not.  This  made  no  difference  with  two 
satellites  in  the  cluster,  but  it  made  a  large  difference  with  10  satellites  by  further 
decreasing  the  true  error.  With  the  above  changes,  several  of  the  cases  presented  in 
the  Ward  thesis  were  validated  with  single  performance  runs  containing  perfect  data 
and  were  then  followed  by  a  full  Monte  Carlo  analysis.  The  set-up  for  the  Monte 
Carlo  analysis  will  be  presented  in  more  detail  later  in  this  section. 

The  first  test  of  the  filter  was  a  stability  analysis.  The  cluster  contained  two 
satellites,  with  the  dynamics  driving  noise  matrix  Qd  entries  equal  to  zero.  The 
resulting  graph,  Figure  5,  shows  a  well  behaved  filter  with  the  filter’s  estimated 
error  decreasing  and  the  initial  true  error  equal  to  zero.  The  true  error,  as  expected, 
began  to  diverge  once  it  had  reached  a  near  zero  value  because  the  filter  assumed  it 
had  reached  perfection  and  ignored  the  new  measurement  data  coming  in.  This  was 
a  result  of  an  inappropriate  choice  of  Qd  =  0. 

The  next  test  was  to  add  driving  noise  to  prevent  the  filter  from  locking  onto 
a  state  and  ignoring  the  measurement  data.  The  initial  conditions  used  above  were 
rerun  to  match  the  Ward  thesis  with  the  Qd  matrix  entries  set  equal  to  1  x  10~ia  km 
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Figure  5.  Stability  Analysis 


for  position  and  1  x  10-17  ^  for  velocity.  Figure  6,  shows  the  true  error  is  still 
diverging  but  at  a  slower  rate,  as  expected.  The  values  chosen  for  the  Qd  matrix 
were  apparently  not  large  enough. 

The  addition  of  more  satellites  to  the  cluster  was  the  next  test  to  determine 
if  additional  measurements  would  improve  the  true  error.  The  cluster  size  was  in¬ 
creased  to  five  satellites  with  all  other  parameters  remaining  constant.  This  again 
improved  the  true  error,  but  unlike  the  results  of  the  Ward  thesis,  the  estimated 
error  in  Figure  7  did  not  not  lock  onto  a  state  and  converge  to  zero. 

The  final  test  was  to  increase  the  cluster  to  the  maximum  of  10  satellites, 
with  the  initial  covariance  matrix  entries  set  equal  to  their  steady  state  values  of 
1  x  10-6  km  for  position  and  1  x  10-12  ^  for  velocity.  Since  the  filter  seemed  to 
be  tuning  with  the  addition  of  less  driving  noise,  the  entries  of  the  Qd  matrix  were 
decrease  from  2  x  10-9  km  to  1  x  10-1°  km  for  position  and  from  2  x  10-14  ^  to 
1  x  10-15  ^  for  velocity. 

Figures  8  and  9  show  the  true  error  is  no  longer  diverging,  but  the  filter  is  still 
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Figure  8.  Satellite  One  with  Steady-State  Covariance  and  Driving  Noise 
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Figure  9.  Satellite  Ten  with  Steady-State  Covariance  and  Driving  Noise 
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overestimating  its  own  error.  The  graphs  also  show  the  error  for  satellite  1  differs 
from  that  of  the  other  satellites.  This  was  also  discovered  in  the  Ward  thesis.  The 
apparent  explanation  lies  in  the  difference  between  the  number  of  components  in  the 
state  vector  for  satellite  1  and  the  other  members  of  the  cluster. 

With  the  code  running  properly,  the  program  was  transferred  to  another  com¬ 
puter  in  order  to  speed  computation  time  for  the  Monte  Carlo  simulation.  The  code 
was  put  on  an  Elxsi  Superframe  with  an  EMBOS  operating  system  and  running  an 
EMS  VMS  operating  system  emulator.  The  code  was  changed  in  several  areas  to 
account  for  system  architectural  differences.  One  obvious  difference  in  the  machine 
architecture  was  the  number  of  digits  carried  by  the  machine  for  single  precision 
which  is  used  in  the  calculations  of  the  estimated  states.  The  VAX  carried  twelve 
digits  and  the  Elxsi  carries  only  nine.  Since  the  error  between  the  true  state  and 
the  estimated  state  was  approximately  .1  meters,  this  difference  did  not  effect  the 
results.  When  the  code  was  tested  on  the  new  computer,  the  results  matched  those 
previously  obtained  with  a  the  running  time  reduced  from  30  minutes  for  a  case  with 
ten  satellites  to  six  minutes. 

The  next  step  was  to  test  the  filter  in  a  full  Monte  Carlo  simulation.  The 
simulation  was  conducted  by  varying  the  seed  values  used  by  the  random  number 
generator  to  determine  the  placement  of  satellites  wound  the  reference  point.  Ten 
different  seed  values  were  used  and  runs  were  completed  with  each  of  these  seed 
values.  The  outputs  were  then  averaged  together  to  obtain  the  final  results.  The 
average  true  error  and  average  estimated  error  were  computed  as  follows: 


average  true  error  at  each  time  interval  =  y] 

k=l 


£&(*(*)-««(*)) 


where  k  represents  each  element  in  the  state  vector  and  »  represents  the  value  com¬ 
puted  for  each  element  at  that  time  interval,  and 
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10 

average  variance  at  each  time  interval  =  ^ 

i=i 


yja(xi)2  +  cT(yi)2  +  a(zi)2 

loTo 


(65) 


where  again  the  a2V)J  values  are  the  eigenvalues  of  the  position  covariance  matrix. 
The  results  of  the  simulation,  Figures  10  and  11,  show  the  same  pattern  as  the  single 
sample  case. 
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Figure  10.  Satellite  One’s  Averaged  Performance 


The  next  step  in  the  testing  of  the  filter  was  the  addition  of  noisy  data  to  the 
system.  The  data  error  standard  deviation  was  initially  set  at  .01  meters  and  the 
driving  noise  entries  were  set  equal  to  2  x  10"9  km  for  position  and  2  x  10~14  ^  for 
velocity. 

Figure  12  shows  there  is  too  much  driving  noise  in  the  system  which  is  causing 
the  filter’s  computed  error  to  diverge  from  the  true  error.  Through  additional  testing 
the  estimated  error  converged  when  the  driving  noise  entries  were  decreased  to  1  x 
10-12  km  for  position  and  1  x  10-17  —  for  velocity.  The  amount  of  noise  in  the 


28 


ure  12.  Satellite  One  of  Ten  with  .01  Meter  Noisy  Data 


data  was  also  increased  to  a  standard  deviation  of  one  meter  to  insure  that  the  data 
covariance  in  single  precision  was  not  zero  km.  The  resulting  graphs,  Figures  13  and 
14,  show  the  estimated  error  is  no  longer  diverging,  but  the  true  error  is  still  too 
large. 
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Figure  13.  Satellite  One  of  Ten  with  1  Meter  Noisy  Data 


Another  method  of  decreasing  the  true  error  is  by  updating  the  filter  at  smaller 
time  intervals.  The  update  cycle  was  decreased  from  once  every  five  minutes  to  once 
every  three  minutes.  As  expected,  Figures  15  and  16  show  that  the  true  error  has 
decreased  to  a  more  acceptable  level,  but  the  filter  is  still  poorly  tuned. 

Before  a  Monte  Carlo  simulation  was  run  for  the  case  with  noisy  data,  the 
inability  of  the  estimated  error  to  converge  to  the  true  error  and  the  increased  true 
error  in  the  case  of  satellite  1,  were  investigated.  The  desired  response  of  the  filter 
would  be  for  the  true  error  to  match  the  filter’s  estimate  of  the  error.  The  filter’s 
computed  error  in  the  above  runs  stays  periodic  and  seems  to  be  unable  to  match 
the  value  of  the  true  error  with  the  information  provided.  A  hypothesis  for  the  cause 
of  this  limit  and  the  increased  true  error  for  satellite  1,  is  the  inability  of  the  filter 
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Figure  14. 


Satellite  Ten  of  Ten  with  1  Meter  Noisy  Data 
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Figure  15.  Satellite  One  with  Three  Minute  Updates 
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Figure  16.  Satellite  Ten  with  Three  Minute  Updates 


to  observe  a  necessary  state  from  the  measurements  available.  The  hypothesis  was 
tested  by  using  the  Kalman  filters  resident  on  each  satellite  in  the  cluster  to  feed 
their  measurement  data  back  to  satellite  1  for  use  in  the  update  computation.  In 
this  manner  the  filter  on  satellite  1  will  have  additional  data  on  which  to  base  its 
estimate. 

The  purpose  of  this  chapter  was  to  validate  the  results  achieved  by  Ward. 
The  cases  presented  in  the  Ward  thesis  were  rerun,  but  due  to  a  small  error  in  the 
measurement  update  cycle,  the  results  did  not  match  those  previously  obtained. 
The  results  obtained  in  this  chapter  show  that  the  true  error  has  improved  over 
that  obtained  by  Ward  and  requires  smaller  values  for  the  driving  noise  entries  to 
keep  from  diverging.  The  filter’s  computed  error  has  not  followed  this  improvement 
causing  the  filter  to  be  poorly  tuned.  The  next  step  is  to  investigate  further  this 
observability  or  observability-like  problem. 
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IV.  Additional  Data  and  Observability  Test 


4-1  Additional  Filters 

One  method  of  obtaining  additional  data  is  to  use  the  Kalman  filter  described 
previously  on  every  satellite  in  the  cluster.  Then  at  each  time  interval,  satellite  2 
through  s  would  send  all  of  their  measurement  data  to  satellite  1  for  inclusion  in 
the  computation  of  the  total  estimated  state  and  the  position  covariance.  Satellite  1 
would  determine  which  satellite  was  providing  the  measurements.  The  variable  used 
in  the  calculations  to  determine  which  satellite  was  currently  providing  measurement 
data  as  well  as  a  clock  pulse  to  satellite  1  was  termed  view.  For  view  =  1, 2, 3, . . . ,  s  — 
1  the  linearized  measurement  matrix  becomes 


In  order  to  avoid  duplication  of  data,  the  number  of  measurements  sent  to 
satellite  1  would  decrease  with  each  additional  satellite  present  in  the  cluster.  For 
example,  if  there  were  four  satellites  in  the  cluster,  satellite  1  would  update  the 
state  with  clock  pulses  from  satellite  two,  three,  and  four,  satellite  2  would  provide 
measurement  data  to  satellite  1  from  the  clock  pulses  it  received  from  satellites 
three  and  four,  and  finally  satellite  three  would  provide  measurement  data  from  the 
clock  pulse  it  received  from  satellite  four.  The  maximum  number  of  measurements 
available  for  a  ten  satellite  cluster  is  55. 

The  entire  linearizaton  measurement  matrix  will  also  vary  depending  on  the 
satellite  that  is  providing  the  measurement  data.  For  satellite  1,  recall  the  overall 
matrix  appears  as  follows: 
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Hi  Hi  0  ...  0 

H2  0  Hi  ...  0 

•  •  •  •  • 

•  •  •  •  • 

•  •  •  •  • 

H.-i  0  0  ...  Ht_i 

When  the  satellite  providing  the  measurements  changes,  this  matrix  remains  the 
same  size  coming  into  satellite  1  by  using  zeros  as  place  holders.  For  measurements 
from  satellite  2  the  matrix  would  appear  as  follows: 

‘0  0  0  0  ...  0  ' 

0  Hi  Hi  0  . . .  0 

H  =  0  Hi  0  Hi  ...  0  (68) 

0  .  ... 

.0  H.-i  0  0  ...  H.- a. 

where  for  i  =  2, 3, . . . ,  a  —  1 

=  0  0  0]  (69) 

Hi  =  0  0  0]  (70) 

The  filter  on  satellite  2, 3, 4, ...,  a  would  have  the  same  values  as  satellite  1  for 
the  initial  covariance  matrix  Po,  and  driving  noise  matrix  covariance  Qj.  Each  of 
these  filters  would  run  independently,  but  send  their  measurement  data  to  satellite  1 
besides  using  the  data  to  update  their  own  estimate.  The  state  vectors  on  each 
satellite  would  also  mimic  the  state  vector  for  satellite  1  by  having  only  five  elements 
in  the  state  vector  for  the  number  corresponding  to  the  satellite.  For  example,  the 
filter  on  satellite  2  would  have  the  following  overall  state  vector  components  for 
i  =  3,4,  ,.,s: 
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Xj 


xa 

z2 

Xj 

V2 

.22. 


Xi 

A  yi 


and  Xi  — 


2i 


Xi 


(71) 


The  initial  testing  of  this  new  system  was  the  placement  of  filters  on  satellites 
1  and  2  in  a  three  satellite  cluster.  The  initial  conditions  were  the  same  as  those  in 
Figure  8.  The  resulting  graphs,  Figure  17  and  Figure  18,  show  the  filter’s  estimated 
error  has  increased  and  it  is  still  unable  to  match  the  true  error.  The  additional  data 
had  no  effect  on  the  true  error  for  satellite  1. 
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Figure  17.  Satellite  One  with  Measurements  from  Satellites  One  and  Two 
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The  additional  measurements  provided  to  satellite  1  do  not  directly  update 
satellite  one’s  components;  they  work  by  allow  the  filter  to  determine  where  the 
other  satellites  in  the  cluster  are  to  a  higher  degree  and  thereby  allowing  satellite  1 
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Figure  18.  Satellite  Three  with  Measurements  from  Satellites  One  and  Two 

to  determine  its  position  to  a  higher  degree.  Figure  17  and  Figure  18  demonstrate 
that  this  approach  does  not  seem  to  be  solving  the  problem.  The  next  step  was  to 
look  at  the  filter  on  satellite  2  and  determine  its  performance.  The  filter  on  this 
satellite  demonstrated  the  same  pattern  as  that  seen  with  satellite  1  by  not  being 
able  to  lock  onto  its  own  state  to  the  same  degree  that  it  locks  onto  other  members 
of  the  cluster.  This  would  seem  to  give  credence  to  the  statement  that  the  change  in 
the  true  error  is  due  to  the  difference  in  the  number  of  elements  in  the  state  vector 
between  the  satellite  the  active  filter  is  residing  and  the  other  satellites  in  the  cluster. 

In  order  to  determine  what  information  the  measurements  were  providing  to 
the  filter,  the  value  of  each  state  position  component  error  was  plotted  against  this 
same  error  after  the  last  update  of  that  time  interval  for  a  one  filter  system.  The 
component  error  is  defined  as: 
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(72) 


e  = 


(x  -  it) 

A  yi 

L  (z  ~  zt) 


Figure  19  and  Figure  20  show  the  results  of  this  test  for  satellite  1  in  a  ten 
satellite  cluster.  The  graphs  show  the  value  of  the  error  changes  very  little  after  the 
update  and  in  some  cases  the  error  after  the  update  has  increased. 
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Figure  19.  Estimated  Minus  True  i  Component  Values 
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This  same  test  was  performed  for  satellite  2.  The  resulting  graphs,  Figure  21 
and  Figure  22,  show  the  value  of  the  component  error  is  decreasing  after  the  update 
as  expected. 

All  of  the  other  members  of  the  cluster  showed  results  similar  to  those  of 
satellite  2.  The  value  of  the  gain  was  checked  to  insure  it  was  not  zero  and  thereby 
not  allowing  satellite  1  to  update.  The  opposite  situation  was  found.  The  update 
values  for  the  components  of  satellite  1  are  higher  than  the  update  values  for  satellite 
2.  Since  the  value  of  the  filter’s  internally  computed  state  component  is  initially 
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Figure  22.  Estimated  Minus  True  z  Component  Values  for  Satellite  Two 

greater  than  the  component  value  from  the  truth  model,  the  greater  update  value 
is  actually  increasing  the  error  in  the  component  rather  than  decreasing  it.  The 
value  of  the  update  can  be  lowered  by  decreasing  the  initial  value  of  the  covariance 
P0  for  satellite  1,  which  in  turn  decreases  the  gain.  With  the  initial  entries  for 
the  position  covariance  decreased  from  1  x  10-e  km  to  1  x  10~17  km  the  case  was 
rerun.  The  results  show  the  gain  decreased  slightly  for  one  period  but  then  returned 
to  its  original  pattern.  One  possible  explanation  for  these  results  is  that  smother 
component  in  satellite  one’s  state  vector  is  unobservable. 

4.2  Observability 

The  inability  of  the  estimated  z  and  x  components  of  satellite  1  to  update 
suggests  that  the  z  and  z  components  of  each  satellites  state  vector  may  be  unob¬ 
servable,  along  with  a  component  used  in  the  propagation  of  x  component.  The  first 
test  was  to  remove  the  z,  i  and  y  components  from  satellite  one’s  state  vector  and 
replace  these  components  in  satellites  2  — s  with  their  relative  components.  The  same 
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steps  that  were  used  to  remove  the  y  component  from  the  system  were  repeated.  For 
i  =  2, 3, s 


A  z  =  zi  —  Zi 


(73) 


A  y  =  yi  -  ji 


(74) 


A  i  =  i,  -  (75) 

For  i  —  1,2,3,  ...,s,  the  new  state  vectors  appear  as  follows: 


xx 


Xi 
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and  Xi  » 


A  Zi 

Xi 


A  in 


A  i 


The  $  matrix  for  satellite  1  was  reduced  to  a  2-by-2  matrix: 


(76) 


4  —  3  cos  ip 
3tj  sin  ip  cos  ip 


(77) 


The  $  matrix  for  satellite  t,  where  i  =  2, 3, ...,  s  remains  a  6-by-6  matrix  with 
the  sign  of  the  x  and  x  components  changed  is: 
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(78) 


0  cos  rp 


where  again  V'  is  rjSt,  and  St  is  the  sampling  time.  The  off-diagonal  matrix  $tl  has 
additional  terms  due  to  the  removal  of  y  term. 


0  0 

6(sin  ip— ip)  ■■  ( cos  ip  —  1 ) 

0  0 

0  0 

677(0)8^-1)  -2sin  \p 

0  0 


(79) 


The  overall  $  matrix  has  the  same  form  that  was  presented  previously  with  the 
$1  and  $,1  matrices  reduced  from  5-by-5  and  6-by-5  to  2-by-2  and  6-by-2  matrices 
respectively. 


This  system  of  state  components  was  tested  with  the  same  conditions  as  those 
present  in  Figure  8.  The  resulting  graphs,  Figure  23  and  Figure  24,  show  a  large  in¬ 
crease  in  the  true  error  for  all  of  the  states.  The  filter’s  estimated  error  has  remained 
very  small  and  is  represented  by  the  dark  line  at  the  bottom  of  the  graphs. 

Since  the  state  vector  for  satellite  1  contains  only  the  x  position  component 
and  still  contains  large  values  for  the  true  error,  the  most  likely  cause  of  the  large 
error  was  the  removal  of  the  y  component  from  the  state  vectors.  This  was  shown 
even  more  clearly  when  the  y  component  was  put  back  into  the  propagation  of  x  and 
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just  the  z  and  z  components  of  the  state  vector  were  removed.  The  resulting  graphs 
Figure  25  and  Figure  26,  show  satellite  1  still  has  the  error  in  the  x  component  of 
its  state  vector  and  the  periodic  estimated  error  since  it  is  not  updating  and  being 
affected  by  the  noisy  data.  The  results  for  satellites  2  —  s  are  closer  to  those  obtained 
with  the  z  and  i  components  remaining,  but  their  removed  has  not  improved  the 
filter’s  performance. 
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Figure  25.  Satellite  One  Error  with  z  and  z  Components  Removed 


Several  other  combinations  of  states  revealed  quickly  that  the  y  component  is 
the  only  state  that  is  obviously  unobservable.  Since  the  error  in  the  system  depends 
on  the  position  in  the  orbit,  the  error  defined  as  the  estimated  state  minus  the  true 
state  of  the  system  was  computed  by  propagating  the  state  vectors  back  to  t  =  0 
after  every  propagation  interval.  The  state  is  propagated  backward  in  time  by  using 
the  inverse  of  the  $  matrix  in  the  following  manner: 

j. 

*o(+)  =  5Z$-1x(+)  (80) 

i=i 


43 


Figure  26.  Satellite  Ten  Error  with  z  and  i  Components  Removed 

I 


To  insure  the  error  propagating  backwards  remained  comparable  to  the  forward  error, 
the  state  vectors  were  propagated  backwards  in  same  number  of  steps  as  it  had  been 
propagated  forward  in  time  through  the  use  of  a  counter.  The  results  of  this  test 
proved  interesting.  The  conditions  used  for  the  run  were  a  two-satellite  cluster  with 
the  Qd  matrix  entries  set  equal  to  1  x  10-10  km  for  position  and  1  x  10-15  for 
velocity,  five  minute  update  cycles  and  no  noise  added  to  the  measurement  data. 
The  resulting  graphs,  Figures  27  -  36,  are  the  exact  opposites  of  each  other. 

The  results  of  these  graphs  again  indicate  the  filter  cannot  observe  some  or  all 
of  the  components  of  the  state  vector,  but  it  can  see  linear  combinations.  The  linear 
combination  of  a  two  satellite  cluster  seems  to  be  the  sum  of  the  components  rather 
than  their  difference.  This  is  the  opposite  of  what  was  expected. 

An  observability  test  was  performed  to  determine  if  one  or  more  of  the  eigen¬ 
values  were  small  enough  that  the  associated  states  can  be  considered  unobservable. 
The  observability  test  was  to  calculate  the  observability  Gramian  matrix,  which  has 
a  discrete  time  representation  represented  by  the  following  equation  (3:243): 
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Figure  32.  Satellite  Two  Estimated  Minus  True  x  Component  Values 
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Figure  34.  Satellite  Two  Estimated  Minus  True  y  Component  Values 
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Figure  37.  Satellite  Two  Estimated  Minus  True  Ay  Component  Values 

Md(0,AT)  6  £*T(i,0)HT(i)RflH(i)4>(i,0)  (81) 

«=1 

The  eigenvalues  of  the  Gramian  at  different  points  along  the  orbit  were  computed  to 
determine  if  there  was  a  group  of  small  values  that  would  indicate  unobservability. 
The  results  of  this  test  showed  that  no  eigenvalue  was  absolutely  zero  which  would 
imply  that  the  covariance  of  that  component  in  the  state  vector  was  definitely  ap¬ 
proaching  infinity,  but  the  spread  of  the  eigenvalue  exponents  made  all  but  three 
of  the  values  untrustworthy.  See  Appendix  B  for  the  exact  values.  The  test  was 
rerun  with  the  measurements  and  $  matrix  in  double  precision  to  attempt  to  get 
more  accurate  results.  Even  in  double  precision  where  the  confidence  in  the  val¬ 
ues  obtained  reaches  eight  decimal  places,  the  eigenvalues  obtained  still  exceeded 
this  spread  and  again  the  results  imply  that  only  three  states  can  be  accurately  de¬ 
termined.  The  numerical  stability  of  the  matrix  was  determined  by  computing  its 
condition  number(3:300).  The  condition  number  of  a  matrix  is  defined  as: 
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*(A)  =  2= 

O’min 


(82) 


where  a^ai  and  a^-n  are  the  maximum  and  minimum  eigenvalues  of  ATA.  When  the 
condition  number  for  the  observation  Gramian  matrix  was  computed,  the  maximum 
and  minimum  eigenvalues  were  on  the  order  of  1025  and  10-8,  which  results  in  a 
condition  number  of  1016.  Studies  show  numerical  difficulties  will  be  experienced  as 
the  condition  number  approaches  10^  where  N  is  the  number  of  significant  figures, 
which  for  this  study  was  16.  The  final  conclusion  resulting  from  the  observability  test 
is  that  the  filter  is  ill-conditioned  and  eight  of  the  states  are  essentially  unobservable 
since  they  cannot  be  determined  by  the  filter  and  are  grouped  together  as  very  small 
numbers. 
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V.  Conclusion 


The  U-D  factorization  filter  under  initial  testing  met  the  required  accuracy  in 
all  the  test  cases.  This  initial  testing  however,  did  reveal  the  performance  of  the  filter 
was  different  for  the  satellite  with  the  functioning  filter,  which  was  satellite  1  for  this 
study.  Under  additional  testing,  this  discrepancy  was  discovered  to  be  caused  by 
the  inability  of  the  filter  to  correctly  update  the  components  of  the  state  vector  for 
satellite  1.  Additional  measurements  were  sent  to  the  filter  residing  on  satellite  1 
but  these  measurements  failed  to  provide  satellite  1  with  the  necessary  information 
required  to  update  its  position  components.  Observability  checks  revealed  none  of 
the  states  were  definitely  unobservable,  but  eight  of  the  position  components  were 
small  enough  that  they  could  be  considered  essentially  unobservable.  Even  though 
the  filter  has  the  most  difficult  time  determining  the  position  of  satellite  1,  the 
inability  of  the  estimated  error  to  approach  the  value  of  the  true  error  for  all  of 
the  satellites  reveals  the  filter  cannot  accurately  determine  the  position  of  any  of 
the  satellites  in  the  cluster.  Currently,  the  relative  measurement  data  provided  to 
the  filter  through  the  use  of  the  onboard  clocks  does  not  supply  the  filter  with  the 
necessary  information  to  determine  the  relative  position  from  a  moving  reference 
point  represented  by  Figure  38. 

The  true  error  of  the  filter  is  within  the  required  accuracy  due  to  initial  values 
given  to  the  standard  deviation  of  the  covariance.  This  is  am  artificial  means  of 
keeping  the  true  error  low.  The  filter  is  updating  the  components  it  can  observe 
and  leaving  the  remaining  states  near  their  steady  state  values.  These  steady  state 
values  are  accurate  enough  for  testing  against  a  truth  model  based  on  the  two-body 
equations  of  motion,  but  if  the  filter  was  placed  in  an  actual  space  environment  with 
many  more  disturbance,  the  inability  of  the  filter  to  update  several  of  the  states 
would  cause  the  filter  to  diverge  from  the  true  position  of  the  satellites  beyond  the 
allowable  accuracy. 
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Figure  38.  Satellite  Cluster  Measurement  Data  and  Required  Position  Determina¬ 
tion 

With  the  results  obtained  from  this  study,  it  is  obvious  the  filter  is  not  func¬ 
tioning  properly.  Since  observability  checks  indicate  a  large  number  of  the  states 
are  essentially  unobservable,  the  development  and  arrangement  system  needs  to  be 
re-examined  from  the  beginning.  Two  recommendations  for  improving  the  filter 
include  changing  the  units  on  the  position  and  velocity  components  so  that  their 
exponential  values  are  closer  together  and  redefining  the  z  component  of  satellite 
1.  The  first  recommendation  will  reduce  the  spread  of  the  steady  state  covariance 
values  which  in  turn  may  provide  for  a  more  numerically  stable  system.  The  second 
recommendation  is  to  force  the  z  component  of  satellite  1  to  be  in  the  reference  orbit 
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which  will  force  that  component  to  become  observable.  If  these  adjustments  do  not 
improve  the  results  it  may  also  become  necessary  to  redefine  the  system  to  include 
only  relative  state  components  rather  than  absolute  components. 
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Appendix  A.  Satellite  True  Error  Accuracy  Requirements 

A.l  Introduction 

This  appendix  outlines  the  procedure  used  to  determine  the  accuracy  require¬ 
ments  for  the  Kalman  filter.  These  requirements  are  used  to  determine  if  the  esti¬ 
mated  position  of  the  satellite  is  within  the  required  accuracy  for  a  receiver  on  the 
ground  to  obtain  a  cohesive  radar  image. 

A. 2  Accuracy  Requirements 

As  stated  in  the  introduction  to  this  thesis,  the  accuracy  of  the  filter  has  been 
decreased  from  one  quarter  of  a  radar  wavelength  to  one  tenth  of  a  radar  wavelength 
to  insure  the  image  is  clear.  The  maximum  wavelength  is  a  function  of  the  size  of 
the  antenna. 

According  to  the  problem  description,  the  cluster  consists  of  ten  satellites  at 
an  altitude  of  1000  km.  Assuming  each  satellite  has  a  ground  coverage  radius  of 
approximately  480  km,  the  angle  subtended  at  the  satellite  by  the  ground  is  defined 
as  (8:A-1) 


‘“<5>  =  0  <83> 

This  arrangement  is  shown  in  Figure  39.  Defining  the  angle  0  as  the  3-dB  beamwidth 
and  solving  for  this  angle  yields 

03dB  =  2tan-1(  j^qq)  =  -895  radians  or  51.282°  (84) 

The  gain  of  the  system  using  the  value  for  the  3-dB  beamwidth  in  degrees  is 

(6:86): 
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Figure  39.  Angle  Subtended  at  the  Satellite  by  the  Ground 


G  =  r,(^)  =  11.61  (85) 

where  rj  is  defined  as  the  antenna  efficiency  and  has  an  assumed  value  of  55  percent. 

If  the  conservative  estimate  is  made  that  the  overall  gain  of  the  cluster  is  just 
the  multiple  of  the  individual  gains  then  (8:A-1): 


Gtot  =  10  x  G  =  116.10 
The  radar  wavelength  can  then  be  found  by 


A*  A. 
Gtot  =  7(“^r) 


(86) 


(87) 
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where  A  is  the  area  of  the  satellite  cluster  defined  as  7r( radius  of  the  cluster)2. 
Rearranging  the  equation  and  solving  for  A  yields 


A  =  2xr(J-) 

&TOT 

For  the  study  presented  in  this  thesis,  the  radius  of  the  cluster  was  .5 km, 
which  when  used  in  the  preceding  equation  yields  the  following  values  of  the  required 
accuracy  of  the  system: 


A  =  216.23m 


(89) 


—  =  21.6m 
10 


(90) 
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Appendix  B.  Observation  Matrix  Eigenvalue  Analysis 

B.  1  Introduction 

This  appendix  provides  the  actual  values  obtained  during  the  eigenvalue  anal¬ 
ysis  of  the  observability  Gramian.  The  initial  conditions  for  the  test  case  were  five 
minute  updates,  steady-state  covariance,  no  noisy  data,  perfect  initial  conditions, 
and  driving  noise  entries  equal  to  1  x  10-10fcm  for  position  and  1  x  lO-15^  for 
velocity.  The  filter  was  allowed  to  run  for  three  orbits  before  the  analysis  was  per¬ 
formed  to  insure  that  the  results  were  not  seriously  impacted  by  the  initial  transients 
in  the  system. 

B.2  Eigenvalue  Analysis 

The  eigenvalue  analysis  was  determined  first  for  the  U-D  factored  form  of  the 
filter  used  throughout  this  thesis.  A  second  filter  used  to  validate  the  results  and 
check  the  programming  of  the  first  filter  was  coded  using  the  standard  Kalman  filter 
propagation  covariance  Equation  (91)  and  the  extended  covariance  update  Equation 
(92)  with  the  entire  filter  in  double  precision.  If  the  factorization  is  correct  the  two 
results  should  closely  match  and  validate  the  conclusion  obtained  in  this  thesis. 


P(-)  =  $P(+)$  +  GdQdG*  (91) 

P(+)  =  P(-)  -  KHP(-)  (92) 

where  K  is  defined  as 

K  =  P(-)Ht  [HP(-)Ht  +  Rf] _1  (93) 
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The  covariance  values  from  each  of  the  filters  were  checked  against  each  other  in 
a  two-satellite  cluster  and  they  matched.  With  the  results  of  the  new  filter  verified  the 
observability  Gramian  was  computed.  The  resulting  eigenvalues  for  the  observability 
Gramian  in  a  two  satellite  cluster  running  for  three  orbital  periods  were  as  follows: 


‘  -.246  x  10-9  ' 

‘-.265  x  10'12' 

-.183  x  10-10 

-.522  x  10“14 

-.373  x  10~14 

-.295  x  lO"14 

-.140  x  10-1S 

-.148  x  10"2° 

.209  x  10"16 

-.148  x  lO"20 

.842  x  10"14 

and  Standard  Filter 

.228  x  10-15 

.710  x  10"12 

.689  x  lO'11 

.591  x  lO"10 

.743  x  lO"11 

.903  x  104 

.903  x  104 

.414  x  107 

.414  x  107 

.717  x  107 

.717  x  107 

(94) 


The  test  run  for  the  first  filter  indicates  that  the  first  eight  eigenvalues  are 
outside  the  eight  decimal  place  limit  from  the  maximum  exponential  value  that  cam 
be  trusted  with  double  precision  accuracy.  The  last  three  values  are  the  only  values 
that  are  within  this  limit.  These  results  were  validated  with  the  second  filter.  The 
agreement  and  accuracy  of  the  last  three  eigenvalues  indicates  that  these  are  the 
only  states  that  the  filter  can  accurately  represent.  Since  the  first  eight  eigenvalues 
are  very  small  they  can  be  considered  essentially  zero  and  represent  the  (essentially) 
unobservable  states. 
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